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Pairing between spinless fermions can generate Majorana fermion excitations. Such excitations 
may exhibit intriguing properties arising from non-local entanglement, including anyonic braid statis- 
tics, teleportation, and enough stability to encode quantum information. But simple models indicate 
that non-local entanglement between Majorana fermions becomes unstable at non-zero tempera- 
tures. We address this issue here by showing that anisotropic interactions between dipolar fermions 
in optical lattices can be used to form domains that significantly enhance thermal stability. We 
construct a model of oriented dipolar fermions in a square optical lattice. We explicitly compute 
the correlation functions defining entanglement. We find that domains established by strong inter- 
actions exhibit enhanced entanglement between Majorana fermions over large distances and long 
times even at finite temperatures. Our approach can be generalized to a variety of configurations 
and other systems, such as quantum wire arrays. 



Optical lattices offer unprecedented tunability in ma- 
nipulating quantum degenerate gases into complex quan- 
tum states [U |2]. The wide variety of lattice geome- 
tries accessible even allow enough control over anisotropy 
to explore dimensional crossover, for example between 
square and chain lattices. But recent developments in 
the cooling of molecules [3-8 and magnetic atoms [9] 
imply that a further level of anisotropy, anisotropy in 
inter-particle interactions, will be tunable in the near fu- 
ture. Such particles have dipolar interactions that, in 
combination with the wealth of lattice geometries possi- 
ble, imply the ability to explore some of the most elu- 
sive yet compelling quantum states, entangled Majorana 
fermions. 

Seminal lattice models demonstrate particle-like exci- 
tations that behave as Majorana fermions thanks to pe- 
culiar non-local symmetries [lOiiUj. A Majorana fermion 
is its own antiparticle. But as excitations in many-body 
lattice models, they exploit symmetries with respect to 
operator strings that entangle with each other over large 
distances to signal underlying topological order, order 
that only responds to non-local perturbations such as the 
surface topology. Pairs of Majorana fermions entangled 
with string correlations support fascinating properties: 
The robustness of these strings motivated proposals for 
topologically protected qubits ^11^ JJ . The crossing of 
string operators is responsible for unusual anyonic braid 
statistics [TTJ[T3]. And string operators connecting these 
excitations also underly theories of quantum state tele- 
portation [T4l[T5lj. 

The zero temperature properties of models hosting 
topological order set the stage for work connected to ex- 
periments. Kitaev's two dimensional toric code Hamil- 
tonian [TT] motivated early proposals in optical lattices 
containing ultracold atoms \i6\ , polar molecules [T7| , and 
atoms in Rydberg states [I8l [TOj . But the one dimen- 
sional Kitaev chain model [10] is one of the simplest mod- 
els supporting free Majorana fermion excitations. Antic- 




FIG. 1: Schematic of dipolar fermions (spheres) in the corru- 
gated potential energy landscape defined by a two dimensional 
optical lattice. The arrows on each sphere indicate that the 
dipoles orient parallel to an applied field, at an angle with 
respect to the x-axis (grey arrow) . The optical lattice imposes 
anisotropic hopping when the lattice depth is large along the 
y direction. The inter-fermion dipolar interaction is repulsive 
along the y direction, but attractive along the x direction for 
small 6. 

ipation of non-local Majorana fermion properties in one 
dimension led to experimental proposals in both optical 
lattices [2QH2^ and solids [TQ1I23H26]. For example, pro- 
posals to observe Majorana fermions in quantum wires 
hosting a topological superconducting state lead to re- 
cent experiments that have shown promising results that 
can be interpreted as localized Majorana fermion excita- 
tions [27] . But prospects for observing interesting effects 
due to non-local entanglement of Majorana fermion pairs 
over long times and large length scales hinge on the sta- 
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bility of non-local string correlators \TT that have not 
yet been directly observed. 

String correlation functions in simple Major ana 
fermion lattice models are unstable at non-zero temper- 
atures. For example, direct calculations show that, in 
spite of a spectral gap, the string correlation functions 
supporting entangled Majorana fermions in Kitaev's two- 
dimensional toric code model vanish at long times and 
large distances because of thermally generated excita- 
tions [TSl [28]. Recent work also argues that Majorana 
fermions in lattice models of topological p-wave super- 
conductors are much more sensitive to thermal fluctua- 
tions than previously thought [29l[30]. A general theorem 
[28^ shows that models of topological order must satisfy 
strict criteria for long-range entanglement to remain re- 
silient against thermal fluctuations at long times. Fortu- 
nately, recent calculations [STJ [32 show that proximity 
coupling [33^ to a reservoir significantly helps in stabiliz- 
ing topological superconductors in wires against thermal 
fluctuations [23ti26]. 

We propose that spatially anisotropic dipolar interac- 
tions in optical lattices [34] offer a powerful tool to ad- 
dress the problem of thermodynamic stability of string 
correlators in lattice models of Majorana fermions. We 
argue that anisotropy in both the lattice and dipolar in- 
teractions can be used to electrostatically copy string 
correlations thus forcing excitations to form arrays of 
strings which we call domains in this work. We solve 
a simple model of dipolar atoms in an anisotropic square 
optical lattice with numerically exact quantum Monte 
Carlo (QMC) simulations to demonstrate that domain 
formation in electrostatically coupled Kitaev chains sig- 
nificantly enhances the built-in entanglement of string 
correlators. Our approach offers a powerful method to 
protect the long sought properties of Majorana fermions 
against thermal fluctuations at the level of the system 
to thereby complement proposals to stabilize Majorana 
fermions with reservoir coupling. 

Model: We first consider a Hubbard model of dipolar 
fermions interacting in an optical lattice and then nar- 
row our analysis to a specific parameter regime. Dipolar 
interactions can arise from an electric dipole moment in 
a polar molecule (e.g., ^^K^^Rb [5 ) or an intrinsic mag- 
netic dipole moment in a single atom (e.g., ^^^Dy [9 ). 
We assume that the dipoles orient along an applied exter- 
nal field (Fig. [T]) but are allowed to hop between nearest 
neighbor (NN) sites. We also assume that a large optical 
lattice depth along the y direction strongly suppresses 
hopping in the y direction. We construct a Hubbard 
model capturing the essential features of this scenario: 

+ [^x{0)nijni^ij + VyfiijUij^i - /J^oTiij] ,(1) 



where we use i = 1, 2, • • • , L^^ — 1 and j = 1, 2, • • • , to 
denote summation indices in the x and y directions, re- 
spectively, a] J creates a spinless fermion at the site (z, j) 

and Hi J = j^i j is the particle number operator. tx{ty) 
is the hopping energy between NN sites in the x{y) di- 
rection. The chemical potential, /io, controls the average 
density. Our approach will work in many different lattice 
geometries of different dimensions but to demonstrate an 
example, electrostatically coupled Kitaev chains, we con- 
sider the square lattice with L = = Ly. (The sup- 
plementary information section discusses boundary con- 
ditions used in calculations.) 

The interaction terms capture the short-range Hub- 
bard part of the dipolar interaction. We will focus on lim- 
its with an energy gap and can therefore truncate the in- 
teraction to NN sites provided that renormalization from 
the long-range part of the interaction yields corrections 
that are much smaller than the gap. We approximate the 
Hubbard parameters with their bare values. The dipolar 
interaction between NN sites in the x and y direction is 
given by F^(l9) = D'^(l - 3cos2 i9)/rg and Vy = D'^/r^ 
respectively. Here tq = A/2 is the lattice constant de- 
fined by the optical lattice laser wavelength, A. In the 
case of electric dipoles, D'^ = p^/47reo with eo the dielec- 
tric constant and p the electric dipole moment. Rotation 
of an applied external field with respect to the optical 
lattice plane tunes the angle that the dipoles make with 
the X-axis, 0. We consider angles such that the fermions 
are attractive along the x-rows (rows along the x direc- 
tion in the optical lattice), 14 < 0, but repulsive along 
^/-columns (columns along the y direction in the lattice), 

Vy>0. 

For a range of yielding Vx < the ground state 
of Eq. ([T]) exhibits p-wave pairing. For tx = ty func- 
tional renormalization group [35] and mean- field [36l [37] 
calculations show a BCS paired state for dipolar interac- 
tions consistent with Eq. ([T]) [29 . p-wave pairing between 
neighbors along x-rows can be modeled by a real space 
pairing field: exp{i^ij)\A\{al^^jalj^a.ja.^^j). But in 
the ty <C tx limit the system can be thought of as weakly 
coupled Luttinger liquids in x-rows. A Luttinger liquid 
theory analysis shows that weakly coupled ID dipolar 
systems also posses p-wave paired order with algebraic 
decaying pairing correlations |38] . Increasing ty connects 
the coupled- ID and 2D square lattice limits. 

ty establishes a Josephson coupling between x-rows. In 
the limit that the inter- x-row tunneling is much smaller 
than the pairing gap within an x-row, we expect second 
order tunneling between paired states: ~ —ty cos{^ij — 
Because the system is otherwise insensitive to 
changes in the relative phase between chains, the Joseph- 
son coupling will tend to align the phase of the pairing 
field between each x-row, j — 0. In what 

follows we assume a uniform pairing field to motivate a 
model that allows the identification of domains of Majo- 
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rana fermion strings with significantly enhanced thermal 
stability. We then discuss corrections in the pair field 
due to fluctuations. 

Mean-Field Limit: We perform a mean-field decoupling 
of Eq. ([T]) to establish a more general model. We take 
ty to be energetically negligible. By decoupling just the 
attractive dipolar interaction term in the x direction we 
find: 



At the Hartree-Fock level the chemical potential is renor- 
malized to /i = /io + '^{'ni,j)\Vx\ — \Vy and the hop- 
ping becomes t = t^ - \Vx{0)\{al^^ja-j). In Eq. (2|, 
we, without loss of generality, tuned to match the 
renormalized pairing term with the hopping by setting 
+ Major ana fermions 

can arise away from this particular point in parameter 
space but this choice eliminates the fermion sign prob- 
lem in our model. This choice also implicitly includes 
the second order effect of the ty term by setting j — 
in the pairing field, t is our energy unit. 

Eq. Q establishes a parameterization of the dipolar 
problem that forms the centerpiece of our work. We will 
use both mean field theory (MFT) and QMC on this 
model to explicitly demonstrate the robustness of do- 
mains against thermal fluctuations. We choose /i = 
(unless stated otherwise) to work near fillings with an 
average of one particle for every two sites. At this filling 
our findings also apply ioVy<^^ since both cases are 
equivalent under a particle-hole transformation on one 
sublattice. 

Eq. ([2]) reduces to known models in certain limits. For 
tl^y — the model does not support kinetics. The 
ground state is essentially a classical state of indepen- 
dent one-dimensional charge density waves aligned along 
^/-columns. In the opposite limit, tjVy 00, we find 
a series of one dimensional Kitaev chains oriented along 
x-rows. For Vy = these x-rows can be mapped onto a 
model of free Majorana fermions. It is known that in- 
teractions within one or two chains can increase the gap 
[39, 40 . But by coupling these chains with strong in- 
teractions, Vy > At^ we will show below that Majorana 
fermions form in unison along columns, to offer a quali- 
tatively more robust mechanism for enhancing entangle- 
ment between entire ^/-columns of Majorana fermions. 
Mapping to Majorana Fermions: We can transform the 
above interacting spinless fermion Hamiltonian into an 
interacting Majorana fermion model by introducing two 
Majorana fermion operators, C2i,j and C2i-i,j, for each 
site of the lattice, (z, j) [10 . These operators satisfy the 
usual fermion anticommutation relations: ^c^, ^, = 



— c^/ ^/C^ ^, for {a, 7^ p'}. But these particles are 

also their own antiparticle: ^ =4/3 (^a z?)^ ~ 
The absence of kinetics along the y direction implies that 
each particle can be labeled with a specific x-iow index, 
j. The Majorana fermion operators then relate to the 
physical fermion operators by a complex superposition: 



{C2i 



ic2i,,)/2 



(C2i-l,i - iC2ij)/2- 



(3) 



We can now demonstrate the existence of edge states by 
mapping Eq. ([2| to Majorana fermion space: 



Hm = 'it ^ C2i,jC2i+l,j + — ^ C2i-ljC2i 



— ^ ^ C2i-l,jC2i,jC2i-l,j + lC2i,j+l. (4) 

Here we see that the first two terms equate to the Kitaev 
chain [the first two terms in Eq. ([2])] and define a free 
Majorana fermion theory. States defined by the dangling 
operators, cij and C2L^,j, at the ends of each x-iow es- 
tablish two-fold degenerate Majorana fermion states that 
can be entangled at T = 0. By non-local entanglement 
we mean that this operator pair hosts a single real dipole 
in a superposition of states that stretches across the j^^ 

X-TOW. 

The Hilbert space of Eq. ([2| also maps onto a quan- 
tum compass model which possesses a spectral gap, AE, 
above the degenerate states for the parameters we con- 
sider here [41, 42 . But the entropy gain in the free energy 
cost to create excitations, AE — ksTS^ can overwhelm 
the energy gap in low dimensions. Here S is the entropy 
and /c^ is Boltzmann's constant. 

The thermal stability of entanglement among excita- 
tions depends on the ratio Vy/t. For weak interactions, 
Vy < At^ the system can be thought of as nearly indepen- 
dent one-dimensional x-rows of free Majorana fermions. 
Non-zero T allows the accumulation of entropy through 
the proliferation of Majorana fermion pairs to destroy 
entanglement at long times and large distances (since 
AE ~ (9(1) and S ~ InL^^) even though fluctuations in 
the pairing field are ignored. But we argue that strong in- 
teraction between x-rows, Vy > At^ significantly improves 
the thermal stability of entanglement between edge states 
because they extend along ^/-columns. Strong interac- 
tions therefore require the creation of entire domains 
(with ~ Ly and 5* ~ InL^^) to destroy entanglement. 
[The domain picture is exact in a mean field limit that 
is shown to be an excellent approximation to Eq. Q in 
the supplementary information section.] We thus propose 
that the interesting entanglement properties between y- 
columns of Majorana fermions are much more thermody- 
namically stable than those between pairs of individual 
Majorana fermions. 
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FIG. 2: The thermal expectation value of the string operator 
plotted as a function of an applied global field [defined in 
Eq. for several different system sizes, L, at = 4.8t and 
/i = 0, computed using QMC. The string operator can only 
be ±1 at T = but may vanish for T > 0. The top (bottom) 
panel shows data for a characteristic low (high) temperature. 
In the top panel we see that very small values of h tend to 
polarize each x-row of the two dimensional lattice into the 
same parity sector, +1. The data collapse and very small 
h values suggest that the parity symmetry of each x-row is 
spontaneously broken. The bottom panel shows that at high 
T, large fields are required to occupy the +1 sector. 



Thermal Stability of Majorana Fermion Entanglement: 
The degree of entanglement between edge states prepared 
at z = 1 and i = Lx is captured by a set of Ly string 
operators that stretch across each x-iow: 



(5) 



i=l 



where j = 1,2, ••• , along the y direction. The sec- 
ond equality shows that this operator is equivalent to 
the fermion parity for each x-row. Each operator Pj has 
two eigenvalues, ±1. Pj commutes with Eq. (|2| to estab- 
lish a 2^y degeneracy [41] |42] on our geometry of choice, 
the cylinder. 

The expectation value of the string operators, Pj, act 
as order parameters. Unique values, (Pj) = ±1, can be 
used to define each sector and therefore indicate stabil- 
ity in the non-local entanglement of Majorana fermions. 
But (P) = indicates that thermal excitations destroy 
any distinction between sectors. We compute (Pj) ex- 
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FIG. 3: Top: The susceptibility of the string- string correla- 
tion function, Eq. ([7|), computed using QMC on Eq. ^ for 
several different system sizes at Vy = 4.8t and /x = 0. Tc is 
defined as the location of the susceptibility peak. The string 
operators tend to order along the y direction for T < Tc. The 
inset shows a schematic of an ordered domain with free Ma- 
jorana fermions forming columns at the ends (dashed lines). 
The domains shrink for T > Tc. Bottom: Tc extrapolated 
to L ^ oo. The solid line is a straight line chi-squared fit 
indicating non-zero Tc in the thermodynamic limit. 



plicitly to show spontaneous breaking of these discrete 
symmetries for Vy > At even at non-zero temperatures. 
To detect the spontaneous breaking of the string symme- 
tries we perturb the above spinless fermion model with a 
weak global field: 



(6) 



The global field, P = L-^ E, =iPj^ imposes a splitting 
between the otherwise degenerate states. We define h = 
hLx with the factor Lx to ensure that the perturbing 
term imposes a non-zero energy splitting per particle, 
between degenerate sectors even in the limit Lx = Ly ^ 
oo. h > favors uniform Pj, (P) = 1. 

We first compute (P) in the limit Vy < 4t using QMC. 
The expectation value is computed using the quantum 
Wang-Landau algorithm in the stochastic series expan- 
sion representation [43, ,44^ (this unbiased method is nu- 
merically exact and is discussed in the supplementary 
information section). For Vy = O.OAt we find (P) ~ to 
within numerical accuracy for h < 0.02t, T = 0.08t, and 
P = 4 — 8 . For Vy = 3.2t we find (P) with increasing 
L. This indicates that the operators defining the sectors 
in one dimensional x-iows alone are extremely sensitive 
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to thermal fluctuations even for small system sizes, as ex- 
pected from the entropy argument above. We note that 
our calculations are time independent. It is possible that 
one can flnd > at short times for weak interac- 

tions. 

We now turn to calculations of (P) in the strongly 
interacting case, Vy = 4.8t, where we expect arrays of 
strings to form stable domains. Fig. [2] shows (P) at low 
and high temperatures computed using Eq. ([6|. Fig. 11 
presents the central result of our work. At high T the 
bottom panel shows that a large value of h is needed to 
stabilize the string operator. But at low T (top panel) 
we flnd that very small flelds tend to force all x-iows to 
spontaneously occupy the lowest energy state in the limit 
h ^ 0. (We do not flnd this behavior for Vy < At.) The 
top panel indicates that ^-columns of Majorana fermions 
located at the ends of the lattice, i = 1 and i = L^^ can be 
prepared in a long lasting entangled state, that stretches 
over large distances, even at flnite temperatures. 

Thermal Stability of Domains: The arrays of string op- 
erators deflning domains are stable at low temperatures 
but eventually break up at large T. To flnd the critical 
temperature at which domain formation becomes unfa- 
vorable, we deflne a string-string order parameter that 
captures the strength of the ordering along the y direc- 
tion: 

The operator O is similar to the static structure factor, 
Sky cx exp [-^ky{j - j')\{njnj'), but with the re- 

placement rijifij' PjPj' and with wavevector ky = 0. 
Thus, the ordering of the string operators along the y 
direction is picked out by the operator O in a way analo- 
gous to deflnitions of density-ordering using peaks in the 
static structure factor. 

We look for long-range order in the susceptibility of O, 
Xo- A peak in xo versus T indicates the critical tem- 
perature Tc at which the large domain breaks up along 
the y direction. For Vy < At we flnd no peaks in our sim- 
ulations and therefore no domain formation for weakly 
interacting chains, i.e., Tc = 0. 

We observe domain formation in xo for Vy > At. The 
top panel of Fig.|3]shows the susceptibility of O computed 
usmg Eq. (H as a function of temperature for Vy = A.8t. 
Above Tc the ^/-columns of Majorana fermions are no 
longer ordered and the average domain size shrinks. At 
these high temperatures the system can be thought of 
as a thermally disordered phase of ordinary fermions. 
The peak also shifts with system size. The bottom 
panel extracts Tc in the thermodynamic limit, yielding 
Tc = 0.25t. Our results agree with studies showing a 
thermal phase transition in the universality class of the 
two dimensional Ising model |45j . 
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FIG. 4: The main panel plots the energy splitting between 
two sectors defined by Pj = ±1 for all x-rows as a function of 
chemical potential for Eq. ([2]) at T = 0.16t and Vy = 4.8t. For 
weak chemical potentials the degeneracy remains robust but 
an L dependent splitting arises for large chemical potentials. 
Inset (a) shows a weak linear increase in density with increas- 
ing chemical potential. The /i ^ slope is quantitatively 
captured by mean field theory. Inset (b) shows a schematic 
phase diagram established by the lifting of the degeneracy, 
horizontal arrow. The vertical arrow indicates the thermal 
phase transition explored in Fig. |3] 



The robustness of the ground state degeneracy also re- 
veals the stability of the string correlations. We denote 
each ground state energy sector by £^(Pi,P2, • • • ). We 
attempt to perturb this degeneracy with a local field, 
a chemical potential shift, to demonstrate the stability 
of the degeneracy that is enforced by the spectral gap. 
The main panel of Fig. |4] shows the energy splitting per 
particle of two different sectors of the Pj operator com- 
puted with Eq. SE = E{-1, -1, • • • ) - ^(1, 1, • • • 
as a function of chemical potential. The flat portion for 
/i/t <C 1 indicates a robust degeneracy. Above /i ~ 1.5t 
the energy splitting acquires a size dependence consistent 
with behavior expected for chemical potentials above the 
gap. Inset (a) shows that the particle density has weak 
linear dependence for /i/t <C 1 which is also captured 
by the MFT discussed in the supplementary information 
section. Our results are consistent with the formation of 
a thermally robust topological phase at T > 0, shown in 
inset (b) of Fig. [ij 

Detection in Optical Lattices: Domain formation can be 
observed directly in time-of-flight measurements. Noise 
correlations between shots of individual time-of-flight im- 
ages relate to the static structure factor, Sk^ of optical 
lattice states [46 . Peaks in noise correlations have al- 
ready been used to identify uniform states in optical lat- 
tices [47l [48] . In the topological phase with domains we 
anticipate the formation of lines, rather than peaks, in 
noise correlations because the density oscillates along just 
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the y direction of the optical lattice for T < Tc. Obser- 
vations of these lines should therefore allow identification 

of Te. 

Entanglement between Majorana fermions could be 
demonstrated through non-local measures. A trapping 
potential, modeled by a spatially varying /i, will allow the 
topological phase near the center of the trap to remain 
intact. But it will leave a topologically trivial phase near 
the edges [large negative fi in inset (b) of Fig. [I]. A dipole 
added to one domain edge should instantaneously occupy 
a superposition of both edge states. Local spectroscopic 
probes (see, e.g., Refs. [20l[22]) applied at each domain 
edge could be adapted to detect the response of one do- 
main edge when dipoles are added to the opposite edge. 
Analogous teleportation experiments have been proposed 
in the solid state |Ml . But ultracold atoms offer more 
direct methods. Recent experiments using high resolu- 
tion spectroscopy to measure particle number parity [49] 
and string operators |50J could be used to explicitly mea- 
sure the string correlation functions, (Pj), for domains 
near the trap center and therefore directly detect long- 
range entanglement. 

Thermal Fluctuations in Pairing: We connected a model 
of oriented fermionic dipoles placed in optical lattices, 
Eq. Q, to a pairing model, Eq. The pairing model it- 
self demonstrates significantly enhanced stability of Ma- 
jorana fermion entanglement (captured by (P)) via do- 
main formation at T > 0. But our specific implementa- 
tion allows thermal fiuctuations of a different quantity, 
the pairing field, exp(i$^^j)|A|, between x-rows. The 
long range component of the dipolar interaction has been 
found to enhance the thermal stability of p-wave super- 
fiuidity [37] but other mechanisms can be engineered to 
further enhance stability. 

Coherent reservoirs can be constructed to suppress 
fiuctuations of the pairing field within the system. Sev- 
eral schemes for inducing a pairing field have been pro- 
posed to generate Majorana fermion models with, e.g., 
dissipation-induced pairing [21] or coupling to a reser- 
voir of bosons [22. These schemes can be adapted to 
ensure pairing fields in dipolar systems. Another alter- 
native could use the proximity effect of a reservoir dipolar 
superfiuid to stabilize the pairing field in the system dis- 
cussed here [3 3 . The pairing gap in the reservoir would 
favor pair tunneling, for weak single-particle tunneling 
between the reservoir and the system, to strengthen the 
pairing field throughout the system. (The supplementary 
information section discusses a candidate optical lattice 
geometry defining a system and reservoir.) Care must be 
taken in engineering the reservoir because excitations in 
the system can be strongly coupled to excitations in the 
reservoir j3l||32j. 

Discussion: We have constructed a mechanism to en- 
hance Majorana Fermion entanglement in an anisotropic 
lattice model. Our approach generalizes to a vari- 
ety of lattice geometries and even other models with 



Majorana fermions provided they take a similar form: 
Ea + Ea,6 Kni^ ^hcrc defines a model with Ma- 
jorana fermions, Fj^J. creates domains with diagonal in- 
teractions between models a and 6, and V^^^ does not 
commute with Hy^ [28 . This class of Hamiltonians also 
applies to Coulomb-coupling in Majorana fermion mod- 
els of quantum wire arrays or quasi- ID tubes containing 
topological superconductors. 
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SUPPLEMENTARY MATERIAL 
Monte Carlo Method and Simulation Parameters 

We use quantum Monte Carlo (QMC) to solve the 
spinless fermion model, Eq. by first using a Jordan- 
Wigner transformation to map the model onto the quan- 
tum compass model [1 . We then perform simulations 
implemented with the Stochastic Series Expansion (SSE) 
[2 combined with the quantum Wang-Landau algorithm 
In this approach the partition function is expanded 
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as a series in powers of /3 = (kBT) ^: 

^ma x /OnTi / ZT\n ^ma x 

Z = Tre-^" = Y: ' ^ = E (8) 

n=0 * n=0 

where A^max is the maximum cutoff of the expansion or- 
der. A/'max determines the lowest temperature that can 
be reached in the simulation, and g{n) corresponds to 
the density of states in a classical system. The distribu- 
tion of g{n) is obtained from a random sampling protocol 
[3| and can be used to estimate the free energy, internal 
energy, entropy, heat capacity, and other properties of 
the system. Note, to measure other physical quantities 
such as the density, density-density correlation, and the 
string operator, we need to accumulate their distribution 
at every order of the series expansion. 

There is a very large energy barrier between different 
sectors defined by the P operator. The large energy bar- 
rier dramatically increases the autocorrelation time in 
conventional QMC simulations with non-local updating. 
Therefore, to reduce the autocorrelation time in QMC 
and enable tunneling between different P sectors, imple- 
mentation of the Wang-Landau algorithm is necessary to 
accurately access the phases discussed in the paper. 

We find that cylindrical boundary conditions (periodic 
in the y direction) are the best choice for this particu- 
lar model. The one dimensional quantum Ising model 
[which can be related to Eq. ^ in the limit Vy = {) af- 
ter a Jordan- Wigner transformation] loses the two-fold 
ground state degeneracy in finite sized systems under 
periodic boundary conditions [H |5] . The degeneracy re- 
mains exact in finite sized systems under open boundary 
conditions [H [5] . The corrections in periodic boundaries 
relate to the fact that the quantum Ising model carries a 
string symmetry. For > in Eq. ([2| we therefore ex- 
pect that periodic boundaries along the x direction will 
lift the degeneracy for finite sized systems. Finite sized 
corrections to the 2^y degeneracy are discussed in Ref. [6] 
for the quantum compass model on periodic boundaries. 
We therefore use open boundaries along the x direction 
(cylindrical boundaries) to prevent perturbations to the 
2^y degeneracy due to finite size effects. 

In the simulation, we check the convergence of var- 
ious physical quantities with respect to the expansion 
cutoff A^max- We find that local quantities such as inter- 
nal energy, average density, density-density correlation 
function, etc., converge much faster than the non-local 
parity operator P at low temperatures, which usually 
requires a much larger A'max- In practice we find the fol- 
lowing values for A^max to be enough for P to converge 
in our simulations in the desired low temperature range, 
A^max = 5000, 8000, and 10000 for L = 4, 6, and 8, re- 
spectively. 



Validating a Mean Field Picture 

To demonstrate the existence of free Major ana 
fermions and domains we perform a mean field decou- 
pling of Eq. ([2| along the y direction. We then verify the 
mean field theory (MET) by direct comparison with an 
unbiased QMC analysis. 

To construct the mean field equations we divide the 
lattice into 2 sublattices, A and 5, along the y direc- 
tion, and decouple the interaction terms. We obtain the 
following 4 coupled mean field equations: 

f^a ^ ^ ^i,a 1 (9) 

where a e {A, B} and [Ia/b = I^a/b - '^yy^'^i^BjA - \)- 
\iA and \iB are chemical potentials for A and B sublat- 
tices, and CLi^j = Fij{alj + a-j)/2. The transformation 
coefficients are given by 

Pi'^i = n 11(1 - 2nk,r) Ha - 2n,^,), (10) 

j'<j k i' <i 

where the summation over j' is extended over all the 
lattice rows except row j. Note that the operator a^^j 
corresponds to a spin ^ operator along the x direction 
in spin space, Sf .^^ based on the Jordan- Wigner trans- 
formation yp. In the spin language, the first equation 
defines a single spin in a magnetic field while the second 
is a quantum Ising model. We use the solutions of both 
of these models [U |5] to solve both models exactly and 
then the coupled equations, Eqs. ([9|, through iteration. 

Solutions of Eqs. ([9| exhibit domains. To see this note 
that the second Hamiltonian in Eqs. ([9|, taken by it- 
self, is a one dimensional model demonstrating Majo- 
rana fermions [7 at zero temperature. It maps onto the 
first two terms in Eq. Q. But the first Hamiltonian in 
Eqs. (|9| just enforces a coupling between x-rows. The in- 
terchain coupling in the original model Eq. ([2| appears as 
a chemical potential shift in the mean field model defined 
by Eqs. 

Eqs. ([9| assume a spatially uniform chemical potential 
(for each sublattice). If this assumption is correct, it 
implies that excitations for any given x-row are copied 
to all other x rows to yield a domain. The existence of 
domains of string operators is therefore implicit in the 
mean field theory but we must validate Eq. ([9| as a good 
approximation to Eq. ([2| to justify this picture. 

We validate Eqs. ([9| by direct comparison with QMC 
solutions to Eq. To compare we compute correla- 

tion functions using both MET and QMC. The following 
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FIG. 5: QMC and MFT comparison of the staggered den- 
sity (top), intra-x-row hopping and pairing correlation func- 
tion (middle), and the inter-x-row density-density correlation 
function (bottom) at Vy = 4.8t. We apply staggered chemical 
potentials /ja and /jb to the A and B sublattices, respectively. 
We choose L = 4 for both QMC and MFT calculations. 



local correlation functions define quantum bond order 
along the x direction and classical bond order along the 
y direction. 



(11) 



Under the spin mapping these correlation functions have 
been studied in a corresponding spin model, the quantum 
compass model [HI . 

Fig. [5] shows that the MFT offers an excellent ap- 
proximation to the QMC results. The large value of 
Vy leads to density- wave ordering along y (large r^). 
But the non-zero values of Tx show quantum correlations 
along the x direction. Therefore both QMC and MFT 
show that the y-columns superpose throughout the lat- 
tice to yield a quantum entangled ground state at non- 
zero temperatures. The good agreement between QMC 
and MFT therefore supports the domain picture implicit 
in Eqs. (IqI). 



System-Reservoir Optical Lattice Geometry 

We show that an optical superlattice can be used to 
host a two dimensional "system" lattice parallel to a two 
dimensional "reservoir" lattice. The system lattice is an 
array of chains in the x — y plane that allow strong tun- 
neling along the x-direction and weak tunneling along 
the ^-direction. The reservoir lattice is a square lattice 
with equal tunneling along both the x and y direction. 



FIG. 6: Plot of the potential defining 
lattice along the z direction for Vz = 
02 37r/2. 



L double well optical 
-15Er, = 0, and 



The increased dimensionality of the reservoir strengthens 
the pair superfluid in the reservoir. A tunable potential 
barrier controls the single particle tunneling between the 
system and the reservoir. 

The optical lattice is formed from three laser beam 
pairs: 1) a double well optical lattice potential, Vzz, 
formed from the interference of counter propagating 
beams along the z direction, 2) a pair of beams with the 
same polarization counter-propagating in the x-z plane, 
to form Vxz , and 3) a similar pair of beams but in the y-z 
plane, to form Vyz. If each beam pair does not interfere 
then the total potential experienced by the particles is: 
Vtotix, y,z) = Vzz{z) + Vxzix.z) + Vyz{y,z). 

The system and reservoir are formed from the double 
well lattice along the z direction. The potential Vzz can 
be formed from the interference of two counter propagat- 
ing lasers with differing wavelengths (See, e.g., Ref p!Q]). 
The distance between the system and the reservoir can 
be changed by using different laser wavelengths to define 
the double well. We choose the wavelengths to differ by 
a factor of 2 to yield: 



[cos {kz — (pi) — cos {kz/2 — ^2)] 



(12) 



Here the wavevector of the primary lattice is k = 27r/A. 
This potential is plotted in Fig. |6] 

We consider an arrangement where the potential es- 
tablished by the remaining beam pairs is given by: 



Vxz{x,z) 

Vyz{y,z) 



Vx [cos ikx) + cos ikz)] 



Vy [cos (ky) + cos (kz)] (13) 
Because the beam pairs forming Vxz and Vyz each have 
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FIG. 7: Plot of the total potential for the system-reservoir 
optical lattice, Vtot- Points are plotted for i;tot < —10^^. The 
parameters are chosen to be: Vz = —15Er, Vx = —0.5Er, 
Vy = -IEr, 01 = -{kn + 27r/1.9), and 02 = -{kn/2 + 
27r/L9). 

the same polarization, they interfere to form a node in the 
z direction at the location of the reservoir. The reservoir 
then experiences a nearly isotropic square lattice even 
with Vx "^y 

Fig. [T] plots an equipotential surface defined by Vtot- 
The potentials are defined in units of the lattice recoil, 



Eji = /2m\^. Here m is the mass of the particles. 
Fig. [T] shows a configuration where the particles in the 
system lattice, near z = 0, have little tunneling along y 
whereas the reservoir lattice, near z = —A, is essentially 
a two-dimensional square lattice. This geometry allows 
a two-dimensional dipolar superfiuid in the reservoir to 
be placed in close proximity to the system lattice. 
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